Mechanistic insights of ABC importer HutCD involved in heme internalization by Vibrio cholerae

Heme internalization by pathogenic bacteria inside a human host to accomplish the requirement of iron for important cellular processes is of paramount importance. Despite this, the mechanism of heme import by the ATP-binding-cassette (ABC) transporter HutCD in Vibrio cholerae remains unexplored. We have performed biochemical studies on ATPase HutD and its mutants, along with molecular modelling, docking and unbiased all-atom MD simulations on lipid-solvated models of permease-ATPase complex HutCD. The results demonstrated mechanisms of ATP binding/hydrolysis and trapped transient and global conformational changes in HutCD, necessary for heme internalization. ATPase HutD forms a dimer, independent of the permease HutC. Each HutD monomer canonically binds ATP in a 1:1 stoichiometry. MD simulations demonstrated that a rotational motion of HutC dimer occurs synchronously with the inter-dimeric D-loop interactions of HutDs. F151 of TM4–TM5 loop of HutC, packs with ATP and Y15 of HutD, initiating ‘cytoplasmic gate opening’ which mimics an ‘outward-facing’ to ‘inward-facing’ conformational switching upon ATP hydrolysis. The simulation on ‘inward-facing’ HutCD culminates to an ‘occluded’ state. The simulation on heme-docked HutCD indicated that the event of heme release occurs in ATP-free ‘inward-facing’ state. Gradual conformational changes of the TM5 helices of HutC towards the ‘occluded’ state facilitate ejection of heme.


Results
Sequence analysis of HutC and HutD. To structurally characterize HutCD of V. cholerae we have performed a BLAST search against PDB with the amino acid sequences of the TMD, HutC (Accession code: C3LWH8) and the NBD, HutD (Accession code: A0A0H3ADP8) and aligned the sequences with their close structural homologs using Clustal Omega (Fig. 1). HutC has shown 45% identity and 67% similarity (for 303aa) with heme transporting permease HmuU of Y. pestis and 45% identity (63% similarity for 339 aa) with BhuU, of B. cenocepacia (Fig. 1a). In addition to the heme transporters, HutC sequence also showed 41% identity and 58% similarity (for 293 aa) with BtuC, the ABC transporter permease involved in vitamin B12 uptake. HutD, the cytosolic ATPase showed a maximum identity of 45% (62% similarity for 242 aa) with HmuV while identity with BhuV of B. cenocepacia was relatively less (36% identity and 51% similarity for 276 aa). Identity of HutD with BtuD of BtuCD was 36% (52% similarity for 217 aa). The significant extent of identity with the heme and B12 ABC transporters of different prokaryotic systems strongly indicates that HutCD is expected to possess a similar overall fold and belongs to the type II class of ABC transporter. The predicted topology of HutC with ten transmembrane helices appears to correspond to the transmembrane segment of the permeases of the type II class (marked in Fig. 1a). Similarly, HutD, like the other NBDs, exhibits sequence and structure conservation across the transporter family 23 and contains several conserved functional motifs such as Walker-A (WA), Walker-B (WB) motifs, an ABC transporter-specific or 'signature' motif, and two shorter sequences containing conserved glutamine and histidine residues (Q-loop and H-loop, respectively) (Fig. 1b).

Cross-linking and size exclusion chromatography suggest dimerization of HutD.
To identify the oligomeric state of HutD, crosslinking experiments were carried out with HutD in the presence and absence of AMP.PNP, the non-hydrolysable analogue of ATP, using Glutaraldehyde as crosslinker (Fig. 2a, Fig. S1). Migration pattern, analysed by 10% SDS-PAGE gel, depicted the presence of dimer (~ 60 kDa) along with monomer of HutD (30 kDa) both in the presence and absence of AMP.PNP (Fig. 2a, Fig. S1). However, tendency of forming dimer is prevalent in the presence of AMP.PNP (Fig. 2a). Little bit of anomaly in movement of proteins in SDS-PAGE (which is observed in case of HutD dimer) is not very uncommon upon cross-linking or chemical modifications of proteins. Moreover. HutD, that was used for cross-linking experiments was pure and functional, and hence the band near ~ 60 kDa has been considered as dimeric band.
To further validate the results of cross-linking experiments, we performed size exclusion chromatography (SEC) using Superdex 200 increase column 10/300 with HutD in free and AMP.PNP treated states (Fig. 2b,c; free state and AMP.PNP treated HutD are shown in red and green colours respectively). In free state, peak II of Fig. 2b having elution volume of 15.36 ml was indicative of monomer along with a trace amount of dimer. However, incubation with AMP.PNP showed a clear shift of peak-II towards dimer with elution volume of 14.68 ml (Fig. 2b,c). This indicates that AMP.PNP binding facilitates dimerization which corroborates with the www.nature.com/scientificreports/ Binding of AMP.PNP with HutD. Participation of A-loop, Walker-A and Walker-B is well known in the ATP binding in the NBDs of type-II ABC importers (PDB codes: 4G1U, 1L2T). Sequence alignment (Fig. 1b) indicated proximity of Y15 of A-loop to the ATP binding site (Fig. 1b). We, therefore, decided to investigate the interactions of AMP.PNP, the non-hydrolysable ATP analog, with HutD by monitoring the intrinsic fluorescence quenching at an excitation wavelength of 280 nm (λ exc = 280 nm, λ em = 295-400 nm) accounting the contributions of tryptophans and tyrosines present at the ligand binding site. Interestingly, with significant fluorescence quenching, HutD has shown considerable interactions with AMP.PNP with a dissociation constant, K d of 94 ± 0.024 μM (Fig. 2d,e). Calculation of stoichiometry indicated binding between HutD and AMP.PNP in a 1:1 molar ratio (Fig. 2f). Negligible quenching upon excitation at wavelength of 295 nm (λ exc = 295 nm, λ em = 308-400 nm), that accounts the contribution of Trp only, ruled out the participation of the Trp residues in ATP binding (Fig. S2) and established that the Trp residues are located beyond the foster distance of ATP binding site of HutD. In turn, these results established the proximity and participation of Y15 in ATP binding.  www.nature.com/scientificreports/ ATP hydrolysis by HutD and calculations of V max and K M . HutD has shown substantial ATP hydrolysis while tested using Malachite green assay (Fig. 2g). 2.5 µM HutD and 100 μM ATP (Sigma Aldrich) were incubated at 298 K for 20 min (Fig. 2g). The release of inorganic phosphate (Pi) was measured at 630 nm upon incubation with Malachite green as per the protocol described earlier 24 . Released Pi from each reaction was quantified by comparing with a Pi standard curve prepared using KH 2 PO 4 (inset of Fig. 2g). HutD without ATP addition served as the control reaction. The experiments were minimally performed in triplicates.
We have determined maximal velocity (V max ) and Michaelis constant (K M ) through time course ATPase assays of HutD. Considering linearity of Pi production during time course experiment, reaction velocity was measured upon incubation of HutD with increasing ATP concentrations (up to 800 μM) through a time scan up to 5 min. Reaction velocity (V 0 ) in terms of Pi release was then plotted against ATP concentrations (Fig. 2h). Fitting of ATPase activities in Michaelis-Menten equation resulted in V max of 3.42 ± 0.10 μM/min and K M of 74.15 ± 9.46 μM for HutD (Fig. 2h).
Identification of pivotal functional residues of HutD through mutagenesis. Based on sequence similarity, we could presume functionally important residues of HutD (Fig. 1b, Table S1). To confirm our propositions, we have prepared four point mutants of HutD namely R18A, K44A, D166A and E167A, and investigated their ATP hydrolysis capacity using Malachite green assays. Time course experiments showed that the amount of Pi produced by wild type (wt) HutD during ATP hydrolysis increased approximately linearly with incubation time for the first 5 min (Fig. 3a). However, the ATPase activity using 100 μM of ATP of the mutants R18A, K44A, D166A and E167A were drastically low compared to wt HutD and remained constant throughout for different incubation times varying from 1 to 5 min (Fig. 3a). The concentration of ATP inside the bacterial host may elevate up to 1 mM under certain conditions 25 . Therefore, we have measured ATPase activities of wt HutD and the variants upon elevating ATP concentrations to 500 μM as well. All four mutants showed negligible ATPase activities even with an elevation of ATP concentration to 500 μM ( Fig. 3a right panel). The increase of incubation time to 20 min also did not affect the ATP hydrolysis pattern of the mutants compared to wt HutD (Fig. 3b).

Modelling of HutD and HutCD in OF and IF conformational states.
Oligomeric states of NBD of ABC transporters and mechanism of ATP binding/hydrolysis differ significantly in the case of type-I and type-II importers. While in the case of the type-I importers dimerization of NBDs are influenced by interactions with permeases, type-II NBDs can form dimer independently 23 . Sequence comparisons proposed HutCD as a type-II importer (Fig. 1). Our biochemical results established that HutD forms stable functional dimer (that can bind and hydrolyse ATP) in the absence of the permease HutC (Fig. 2a-c). Therefore, we decided to investigate the dynamic nature of HutD dimers in the OF and IF states.
The structures of the two type-II Heme importers, HmuUV of Y. pestis (PDB code: 4G1U) and BhuUV from B. cenocepacia (PDB code: 5B57), in two different conformational states, served as templates to prepare the models of HutD and HutCD in the OF and IF states. First, we prepared the models of HutD dimers in OF and IF conformations (Fig. 4a,b). Previous studies on type-II ABC transporters suggested that ATP hydrolysis mainly occurs in the OF state 26 . Therefore, knowledge-based docking of two Mg 2+ -ATP molecules was performed with HutD in OF state (Fig. 4a,c) where coordinates and binding mode of ATP in MJ0796 NBD dimer (PDB code 1L2T) was treated as template. Dimers of HutC in OF and IF states were also modelled in a similar manner. All models were prepared initially using ps2v2 27 and later verified by Alphafold2 28 . HutCD assemblies in IF and OF states were then generated considering the assembly structures of HmuUV and BhuUV (Fig. 4d-i).
Altogether, we have prepared five sets of complexes: (1) Mg 2+ -ATP bound HutD dimer in OF state; (2) HutD dimer in IF state; (3) lipid bilayer embedded Mg 2+ -ATP bound HutCD assembly in OF state; (4) lipid bilayer embedded HutCD in IF state; (5) Heme bound lipid bilayer embedded HutCD in IF state. MD simulations, each of 1 μs, were carried out on all complexes with the Amber-18 and AmberTools-18 software packages 29 . We have used DMPC (1,2-dimyristoyl-sn-glycero-3-phosphocholine) molecules for constructing the bilayer and CHARMM-GUI 30 was used to prepare the solvated protein-embedded lipid bilayer system before simulations. www.nature.com/scientificreports/ An embedded model is shown in Fig. 4i. However, the lipid bilayer and the solvents molecules were not shown in the subsequent figures to maintain clarity.

Dynamics of HutD dimer in free and in HutC bound states.
In order to understand the collective internal motions of HutD dimers in the OF and IF states, we have calculated root-mean-square deviations (RMSDs) of the Cα atoms from the snapshots of the trajectories up to 1 μs for the two chains of NBDs of HutD dimer and HutCD dimer (Fig. 5a,b). In all cases, RMSD values converged well to a constant value after approximately 75 ns. RMSD values of the two HutD protomers (named as HutD1 and HutD2) were plotted separately ( Fig. 5a,b). The RMSD values for HutD1 of HutCD OF, HutCD IF, and HutD IF and OF were steady around 3.0 to 3.5 Å (Fig. 5a). However, HutD2 in OF state showed relatively lower RMSDs ranging around 2 Å (Fig. 5b). This indicates that NBDs acquire stability upon Mg 2+ -ATP binding, although asymmetric movements were observed between the monomers (Fig. 5a,b). B-factors of the two chains of HutD were also plotted separately. The highest B-factors were observed for HutD and the B-factors were lowest for HutD of HutCD in ATP-Mg 2+ bound state (Fig. 5c,d). HutD of HutCD IF state showed intermediate flexibility (Fig. 5c,d). Further analysis of B-factors of the two monomers of HutD  Fig. 1b) are also evident for all the states, which is more distinct in case of HutD1 compared to HutD2 (Fig. 5c,d). As evident from large B-factor values, high flexibility were observed in certain surface exposed regions of HutD (Fig. 5c,d).
Interactions near ATP binding site of HutD dimer. Analysis of the simulation trajectory of HutD dimer in the presence of Mg 2+ -ATP bound OF state revealed that Y15 of the A-loop packs hydrophobically with the adenine base of ATP where R18 interacts with the ribose sugar of ATP (Fig. 5e). During dynamics, Mg 2+ consistently remains bound to E167 of Walker-B and is poised to hydrolyze ATP (Fig. 5e,f). D166, the neighboring residue, interacts with S45 of Walker-A and thereby stabilizes the Walker-A conformation. K44 acts as the pivotal residue of Walker-A that interacts and stabilizes β-and γ-phosphates of ATP (Fig. 5e). These observations corroborate with abrogation of ATP hydrolysis of R18, K44, D166 and E167 to Ala (Fig. 3).
In the OF state, the D-loops of the two HutD monomers interact with each other at the dimeric interface which is mediated primarily by T169 and S170 (Fig. 5f). H200 also participates in dimerization (Fig. 5g). H200 of both the chains interacts with the D-loops belonging to the trans protomer (Fig. 5g). N40 is found to play a www.nature.com/scientificreports/ dual role. This residue is capable of interacting both with the γ-phosphate of bound ATP and D173 of the D-loop belonging to the trans-acting HutD ( Fig. 5f,g). ATP binding was found to be assisted by ' ABC signature motif ' in BtuCDF 26 . However, we have hardly observed any contribution of ' ABC signature motif ' in ATP binding for HutCD ( Fig. 5f).

Dynamics of HutCD in the OF state.
We observed that RMSDs of the two HutD chains of HutCD converged after 75 ns and remained steady around 2.5 Å and 3 Å respectively. To investigate the flexibility and conformational shift of TMDs in HutCD, we have calculated the RMSDs of HutC1 and HutC2 in OF and IF states from MD trajectories run up to 1 μs. RMSD plots of two TMD chains (Fig. 6a,b) indicated that the values converged after approximately 75 ns. While RMSD values of HutCs in the IF state remained stable around 3 Å (Fig. 6a,b), significant fluctuations were observed in the RMSDs of HutCs in the OF state, ranging between 4 and 5 Å. Overall superposition of the snapshots in OF state (ranging between 100 and 750 ns) demonstrated high flexibility of the loops exposed to periplasmic side and encompassing the PLBP binding site (aa 30-55 of HutC) (Fig. 6c). Although these loops of both HutC monomers were apparently in the upright conformation in the staring model (Fig. 4d), after dynamics they culminated in different conformations (Fig. 6c). In one HutC monomer, this loop region was found in upright conformation whereas in the other it acquired an inclined con- www.nature.com/scientificreports/ formation (Fig. 6c). Asymmetric movement of this loop in two HutC monomers might be attributed to ease the accommodation of the PLBP.

Synchronized conformational shift of HutD and HutC in ATP-Mg 2+ bound OF state
HutCD. Analysis of the simulation trajectory of HutCD OF state indicated that the interactions between the D-loops of the two HutD monomers are instrumental in the transmission of conformational shift to HutC. In HutCD, the distance between S170 of the SALD motif of the D-loops of the two monomers varied from 2.9 (103 ns) to 17 Å (606 ns) ( Table 1). Upon proximity, the D-loops interacted with each other through H-bonding between S170 and/or T169 (Fig. 5f). N40 stabilized the γ-phosphate of ATP and interacted with D173 belonging to the 'SALD motif ' of the D-loop of trans-acting protomer (Fig. 5f,g). However, this interaction was inversely related to inter HutD interactions through the D-loops, which is evident from Table 1. Upon proximity of the D-loops of the two monomers of HutD, N40 interacted with the γ-phosphate of ATP practically facilitating ATP hydrolysis. Available structures suggest that a cytoplasmic loop, called L-loop that is located between TM6 and TM7 helices of TMD, shares extensive contacts with the complementary groove on the NBD surface and is proposed to be critical in TMD-NBD communication 31 . Such interactions were evident in HutCD as well (Fig. 4d). Superposition of the snapshots, collected from 50 to 580 ns of the simulation run of HutCD, established retention of proximity and hydrophobic interactions between L-loop of HutC and Q-loop of HutD (Fig. 6d). The conserved Q85 of the Q-loop of HutD was also found to interact with the cis-acting D-loop. The conformational shifts in HutD dimer, caused by the position and the interactions of the D-loops, therefore, relayed to the L-loop of HutC through the Q loop of HutD. Trajectories of HutCD in OF state between 140 and 300 ns showed vigorous structural movements, and after 300 ns there has been a continuous loosening and widening of the overall structure, with concomitant dissociation of HutD1-HutD2 interface by distancing of D-loops ( Table 1).
Opening of the 'cytoplasmic gate'. The 'cytoplasmic gate' that resides in cytoplasmic side of TM5 is known to be instrumental in heme translocation, as evident from the structures of HmuUV and BhuUV 13,14 . F151 that belongs to the loop between TM4 and TM5 of HutC located near cytoplasmic gate is found to be important in terms of transmission of the signal. At the beginning of simulation, this loop was 13 Å away from HutDs. During dynamics, F151 was seen to pack with adenine base of ATP and Y15 of the A-loop of HutD (Fig. 6e). Interestingly, the involvement of F151 in such hydrophobic packing is not consistently simultaneous in each HutC-HutD duo. In early stage of simulation, such as 80 ns to 150 ns, the interaction remained simultaneous in both HutC1-HutD1 and HutC2-HutD2. However, after 200 ns the interactions were not that consistent. Subsequently, F151 of HutC packs with Y15 of HutD but adenine base of ATP remains sufficiently away from those residues. In these stages of the trajectory, such interaction was observed in one HutC-HutD while in the other the TM4-TM5 loop was away from the ATP binding site. At the beginning of MD simulation 'cytoplasmic gate' was closed (Fig. 6f) and as expected with OF structure, periplasmic side of the heme uptake channel was open (Fig. 6g). During the progress of the simulation run, packing of F151 of HutC with ATP and Y15 of HutD increases the distance of TM5 in cytoplasmic side essentially opening the cytoplasmic gate (Fig. 6h). A superpo- Table 1. Random frames taken in order to see any proximities between S170 of HutD1:N40 of HutD2, N40 of HutD1:S170 of HutD2 and S170 of HutD1:S170 HutD2.
Trajectory of HutCD OF state (ns) S170 of HutD1/N40 of HutD2 (Å) N40 of HutD1/S170 of HutD2 (Å) D-loop distance S170-S170 www.nature.com/scientificreports/ sition of the snapshots at 150 ns on that of frame 1 depicts such opening of the 'cytoplasmic gate' (Fig. 6h). This conformational shift made us hypothesise that ATP binding and hydrolysis facilitates opening of the 'cytoplasmic gate' for heme uptake.
Occluded state is observed in the simulation of ATP-free IF state of HutCD. As expected, the HutCD IF model showed inward facing orientation of the TM5 helices (Fig. 7a). Overall RMSD values of HutC and HutD dimers remained around ~ 3 Å during simulation. Throughout the dynamics, the periplasmic side remained tightly sealed by salt bridge interactions between D182 and R186 along with hydrophobic packing between L185 of two H5a helices. Although the crossing angle of H5a-H5a was little lower (50.8°) than the structure of BhuUV (73.1°), like BhuUV, the position of H5a helices remained almost unchanged during dynamics. Notable transition in the orientation of the TM5 helices was observed from 277 ns. At this stage and onward, drastic reduction in the inter TM5-TM5 crossing angle occurred from 23.9° to 6.4° (Fig. 7b). Superposition of the snapshots of HutCD IF, collected at ~ 277 ns, with the occluded ATP bound crystal structure of vitamin B12 importer BtuCD (PDB Code: 4FI3) showed a RMSD of 3 Å (Fig. 7c). Maximum differences were observed among the NBDs, while orientation of the TM5 helices of the two structures were outstandingly similar to each other (Fig. 7c,d). The crossing angle/distance of TM5 axes was 6.4°/16.2 Å in case of HutCD, and that of the ATP bound occluded structure of BtuCD is 11.2°/19.0 Å (PDB Code: 4FI3). These observations suggest that unbiased MD simulation of the HutCD in IF state leads to gradual transition to the ostensible occluded state from ~ 277 ns by exclusive movement of the inter TM5 helices (Fig. 7). The NBD domains of HutCD IF state remained oriented throughout dynamics, where the main flexible portion was the A-loops (Fig. 5c). In IF state, the trans-acting ' ABC signature motif ' was found near the Walker-A, although the distance was not equal in case of the two protomers. Gradual proximilty of Walker-A and www.nature.com/scientificreports/ trans-acting ' ABC signature motif ' with the progress of simulation provided compactness to the NBD dimer. The D-loops at the NBD interface were less flexible than the D-loops of Mg 2+ -ATP bound OF state HutCD. In contrust to the OF state where D-loops of HutDs were found to interact, in the IF state, inter D-loop distance was ranging from 6 to 11 Å in most of the frames.

Identification of heme release pathway in HutCD IF state. Naoe et al. stated that D112 of the TM2-
TM3 loop of the permease BhuU of B. cenocepacia is the only charged residue exposed to the inward facing heme translocation pathway of BhuUV and plays critical role in the heme translocation process, because mutation of this residue to the hydrophobic residues impaired transport activity 14 . Corresponding residue of HutC is E92 and as expected, this residue of both the HutC monomers are exposed to the heme translocation channel in the IF model of HutCD to impart polarity (Figs. 1a, 8a). We, therefore, intended to investigate the binding pattern www.nature.com/scientificreports/ and movement of heme, after it reaches the 'cytoplasmic gate' . For that reason, heme was docked near the cytoplasmic gate of inward facing HutCD using Autodock Vina 32 (Fig. 8b,c). The hydrophobic site was made of L159 and V163 of TM5 along with L90 of TM2-TM3 loop. One propionate group was hydrogen bonded with R87 of the TM2-TM3 loop (Fig. 8d). Soon after beginning of simulation, the heme molecule moved to the pocket between two HutC-HutD dimers ( Fig. 8e-g). At this stage, the propionate group was interacting with the backbone NH groups of the Q-loop residues L89 and T90 of HutD1. Heme stayed inside this pocket 'flip flopping' around for a significant period of time. Subsequently, heme was found to interact with N40 of Walker-A and the D-loop residues of HutD. Polar interactions of the propionates of heme with Q86, S88, T90 of HutD (HutD1) were steady from 350 to 430 ns (Fig. 8f). Additionally, F151 of HutC was found to pack with heme (Fig. 8g). Notably, F151 of HutC was crucial for cytoplasmic gate opening during ATP hydrolysis by HutCD in OF state. The pocket between two HutC-HutD dimers started squeezing around 455 ns. Around 465 ns the heme changed orientation and moved near the other HutD monomer (HutD2) where the propionates were interacting with S138 of HutD2 (Fig. 8g). Heme was on the verge of release to cytosol around 472 ns where heme was in contact with S138 of the ' ABC signature motif ' of HutD2. At 473 ns, the propionate group of heme had polar interaction with S136 followed by H128 and the molecule is finally released to the cytosol (Fig. 8e-i). A detailed list of interactions between HutCD and heme is given as Table S2.
During the process of heme release, no significant structural change occurred in the TMDs except orientation of TM5 (Fig. 8h). After 400 ns, the TM5 helices started acquiring an occluded conformation which was prominent around 460 ns (Fig. 8i). Attainment of occluded state is evident from the superposition of the snapshots collected at 85, 121, 410 and 468 ns (Fig. 8i). Here, conformation of TM5 remained unchanged in HutC1 where cytoplasmic side of TM5 of HutC2 shifted about 5 Å towards central axis of the translocation channel (Fig. 8i).

Discussion
Dimerization of the NBDs in ABC importers of prokaryotes, whether nucleotide driven or independent, is indispensable, since this is a vital early event of ligand translocation. Our results authenticated HutD as a type-II ABC importer, since it forms dimer even in the absence of the permease HutC (Fig. 2a-c). Based on the sequence comparisons, Y15 of the A-loop of HutD was predicted to pack with the adenine base of ATP (Fig. 1b). Substantial fluorescence quenching in the presence of AMP.PNP established significant binding of ATP with HutD in 1:1 stoichiometry with the involvement of Y15 (Fig. 2d-f). While wt HutD was an efficient ATPase, R18A, K44A, D166A and E167A turned out as loss-of-function mutants (Fig. 3a,b). Abrogation of ATPase activities in the mutants together with modelling, docking and MD simulations results suggested that R18 of the A-loop and K44 of Walker-A of HutD stabilize the ribose sugar and the β-phosphate of ATP respectively and hence, ATP binding is compromised upon removal of those side chains (Figs. 1b, 3a,b, 4c, 5e). Based on the sequence analysis (Fig. 1b) we presumed that D166 and E167 are crucial Walker-B residues for conducting ATP hydrolysis event. Our observations demonstrated that E167 is implicated in Mg +2 binding for ATP hydrolysis while neighbouring D166 interacts with S45 of Walker-A to facilitate ATP binding by stabilizing Walker-A. Rather, a network of polar interactions was observed between S45 of Walker-A and D166, E167 of Walker-B as evident from the snapshot of MD simulation (Fig. 5e). Mutation of D166 to Ala presumably generates more space for E167 side chain and at the same time turns E167 more flexible because of impaired polar interactions stated above, eventually disturbing Mg +2 binding/ATP hydrolysis.
MD simulation trajectory of HutD in ATP bound state depicted that D-loops of the two HutD monomers interact at the dimeric interface and H200 critically contributes to the dimerization by interacting with the transacting D-loop (Fig. 5g). MD simulation results further suggested that ATP binding to HutD does not need any assistance from the trans-acting ' ABC signature motif ' (Fig. 5f), like it was observed before in the crystal structure of BtuCD-F 26 . Rather, in HutCD, ' ABC signature motif ' interacts with the heme during its ejection to the cytosol. Presumably, the translocation of bigger ligand vitamin B12 requires a wider channel compared to heme (HmuUV structure) and the NBDs need to dimerize accordingly. The difference in ATP binding mechanisms might be attributed to the differential inter dimeric interactions of BtuD and HutD which are essential for the ligand-and species-specific transportation pathways.
Conformational shifts of certain transmembrane helices, especially the movable helices TM5 and H5a regulate the ligand import processes through ABC importers. Substrate translocations through the translocation channels of type-I ABC transporters MalK of E. coli are dictated by 'tweezer' like or 'twisting' motions of the NBDs 33,34 . For type II importers, substrates translocation by open/closure of the cytoplasmic gate is thought to occur via 'peristaltic' motions of the TMDs 16 . In HutCD, a rotational motion occurs in the TMDs in synchronous manner with the inter NBD D-loop interactions ( Table 1). The more the separation of the inter HutD D-loops, the more rotation of the NBDs which is then transmitted to the TMDs. The distance between the two HutC-HutD dimers was found to be altered accordingly. However, asymmetry has been observed in the interface interactions between HutC and HutD for each HutC-HutD duo.
The interaction of the SKF 151 loop of HutC, located between TM4 and TM5 near the cytoplasmic gate, with ATP and the A-loop (Y 15 GSR) of HutD is one of the most prominent facts that tie the event of ATP hydrolysis with the conformational change in HutCs to facilitate ligand translocation (Fig. 6e). Packing of F151 with the adenine base of ATP and Y15 of the A-loop of HutD leads to the widening of the cytoplasmic gate where the periplasmic side of the TM5 helices remained practically unchanged (Fig. 6e,h). Interestingly, in the later stage of the dynamics (~ 400 ns and beyond), which might be treated as post hydrolysis phase, continuous relaxation of the assembly was observed and, in this situation, TM5 helices prefers to acquire relatively parallel conformation. Inward facing assembly was also culminated to an occluded conformation, similar to that observed before for www.nature.com/scientificreports/ vitamin B12 importer BtuCD (Fig. 7c,d). We, therefore, hypothesize that the occluded state is the most preferred conformation of HutCD. H5a was seen to play a key role as a gating helix to forming the interaction site for the PLBP, BhuT in the case of BhuUV-T assembly 14 . Also, mutation of R176 in H5a of HmuU was reported to decrease affinity to HmuT and decrease heme-transport activity 13 . Although TM5 helices showed visible transition near the cytoplasmic gate, conformation of H5a helices, located at the periplasmic side, practically remianed unchanged during the dynamics. We, therefore, hypothesize that the movement or change in orientation of H5A helices are most likely PLBP dependent, and subsequent separation of the PLBP from the TMDs may result in the transition of orientations of the H5A helices.
The heme release mechanism has been deciphered for the first time which was unknown so far. Our extensive simulation showed that the release of heme is primarily dictated by the interaction of the propionate groups with the polar residues of HutC-HutD dimeric pocket (Fig. 8f,g). Interaction of heme with the residues of D-loop, Walker-A, Q-loop and ' ABC signature motif ' (Fig. 8f,g, Table S2) further suggested that heme release event takes place in the absence of ATP. Because the presence of ATP would have hindered involvement of these motifs in the heme release process. As observed in case of heme-free HutCD in the IF state, heme bound form also acquired occluded state at the verge of heme release. Parallel conformation of TM5 helices in the occluded state probably leads to the squeezing of inter-dimeric pocket which eventually facilitates the heme release (Fig. 9). Hence ABC importers like HutCD prefers to acquire occluded conformation before starting the next cycle of ATP hydrolysis dependent heme translocation (Fig. 9).
Type-II importers HmuUV, BhuUV or BtuCD were predominantly in different states while trapped in the crystal lattice 13,14,26 . Notably, Y15 of HutD or F151 of HutC, which stack with the adenine base of ATP leading to conformational changes in HutCD assembly, are not conserved in HmuV or BhuV (Fig. 1b). Although similar hydrophobic residues are observed in HmuUV in the close vicinity of the respective residues, BhuUV is quite different in this regard (Fig. 1). Furthermore, S170, the crucial D-loop residue of HutD is replaced by Ala in BhuV (Fig. 1b) generating a possibility of compromised D-loop interactions. Collectively, the observations intend us to believe that although there is a common minimum procedure to internalize a specific nutrient, detailed mechanism of uptake is primarily dictated by the requirement of that specific nutrient in the respective species.

Methods
Cloning, site-directed-mutagenesis, overexpression, and purification. Chromosomal DNA of V.
cholerae strain O395 was used as the template to amplify the region encoding Full length HutD of 259 aa (Accession code: A0A0H3ADP8). The 780 bp HutD PCR amplicon and the pET28a + vector (Novagen) with NdeI and BamHI restriction sites were ligated using T4 DNA ligase and the appropriate clones were selected using E. coli XL1-Blue cells with kanamycin resistance. The cDNA of HutD (amino acids 1-259) was PCR amplified and was cloned in pET28a + within NdeI and BamHI restriction sites. The recombinant protein with N-terminal His 6 -tag was overexpressed in BL21 (DE3) by IPTG induction and purified by Ni-NTA affinity chromatography. The construct was verified by restriction digestion analysis and commercial DNA sequencing. The mutants R18A, K44A, D166A and E167A were prepared by two-step PCR amplification using previously cloned WT HutD www.nature.com/scientificreports/ plasmid as template. These mutant amplicons were also cloned in pET28a + vector using NdeI and BamHI as restriction sites. Sequences of the variants were verified by commercial sequencing. These recombinant proteins with N-terminal His 6 tag were overexpressed in E. coli BL21 (DE3) cells in the presence of IPTG.
Overexpression and purification of recombinant proteins were performed as per protocol described by Dey et al. 24 . In brief, 10 ml of LB broth, supplemented with kanamycin, was inoculated by a single colony, and grown by overnight shaking at 310 K. 1 l LB broth with kanamycin was inoculated with 10 ml of the overnight culture and the culture was grown further at 310 K until the OD 600 reached 0.6. The cells were induced by 1 mM IPTG and grown at 310 K for another 3 h. The cells were then harvested at 4500×g for 20 min at 277 K and the pellet was resuspended in 10 ml ice-cold lysis buffer-L (having 50 mM Tris-HCl pH 8.0, 300 mM NaCl, 5 mM MgCl 2 , 10% (v/v) glycerol). PMSF (1 mM) and lysozyme (1 mg/ml) were added to the resuspended solution and it was lysed by sonication on ice. The cell lysate was then centrifuged at (12,000 RPM for 45 min) at 277 K. The collected supernatant was applied onto a nickel-nitrilotriacetic acid (Ni 2+ -NTA) affinity chromatography media (Qiagen) that was previously equilibrated with buffer-L. The 6× His-tagged recombinant proteins were eluted with 50-200 mM imidazole gradient. Subsequently, imidazole was removed from proteins through buffer exchange using Amicon centrifugation units. The wt HutD was concentrated up to 300 μM. All mutants were purified using the buffer containing 50 mM Tris-HCl pH 8.0, 300 mM NaCl, 10% (v/v) glycerol with the same protocol and concentrated up to 250 μM. The homogeneity of the purified proteins was checked using SDS-PAGE with 12% polyacrylamide concentration.

Cross-linking experiments. Crosslinking experiments had been carried out with WT HutD and WT
HutD in the presence of non-hydrolysable ATP analogue AMP-PNP using Glutaraldehyde (Sigma-aldrich) as crosslinker. HutD was purified in cross-linking buffer (50 mM Na-Phosphate pH 8.0 and 150 mM NaCl) followed by incubation with AMP.PNP (20 mM) on ice for 1.5 h. Increasing amounts of glutaraldehyde were added to the mix to a final concentration range of 0.0025-0.04%. The reaction mixture was further incubated at room temperature for 15 min. Crosslinking reactions were quenched by adding 50 mM Tris-HCL pH 8.0, followed by incubation at 97 °C for 10 min in standard SDS loading dye. Reaction products were separated by electrophoresis and analysed on 10% SDS-PAGE polyacrylamide gel.

Gel filtration assay. Size Exclusion chromatography experiments were performed in a Superdex 200
Increase 10/300 GL using AKTA purifier (cytiva). Blue dextran was used to identify the void volume. The column was further calibrated with standard molecular weight calibration kit (GE Healthcare) containing Ferritin (440 kDa), Aldolase (160 kDa), Ovalbumin (45 kDa) and Lysozyme (14.3 kDa). The standard graph was prepared against relative elution volume (V e /V o ) in X-axis [where V e is the elution volume and V o is the void volume] and the log molecular weight in Y-axis. For the AMP.PNP bound state, 250 µM HutD was preincubated for 1.5 h at 4 °C with 500 µM AMP.PNP in the same buffer-L. For analysis of HutD and HutD with AMP.PNP, in each run, 500 μl of protein was injected to the SEC column pre-equilibrated with buffer containing 50 mM Tris-HCl (pH 8.0), 300 mM NaCl, and 5 mM MgCl 2 and ran at a flow rate of 0.5 ml/min. The peak fractions were collected and analysed in 12% SDS-PAGE.
Fluorescence quenching study. As per protocol described by Agarwal et al. 22 , fluorescence measurements were carried out in Hitachi F-7000 spectrofluorometer, using quartz cuvettes of 1 cm path length. Changes in fluorescence of tryptophan and tyrosine residues were measured at an excitation wavelength of 280 nm and the emission spectra were recorded between 295 and 400 nm. Sole contribution of Trp residues upon AMP. PNP binding was also monitored with an excitation at 295 nm and emission between 308 and 400 nm. In both the cases, slit widths was kept consistently at 5 nm. All reactions were carried out at 298 K. The reactions were performed in a buffer containing 50 mM Tris-HCl (pH 8.0) and 300 mM NaCl. Equilibrium titration of HutD was carried out with AMP-PNP as ligand and the changes in fluorescence emission intensity were measured in the presence of increasing concentration of ligand. The concentration of HutD was 2.5 μM and ligand concentrations varied from 0 to 2 mM.
The binding stoichiometry was determined using the protocol described by Mani et al. 35 . The plot of log(F 0 − F)/(F − F ∞ ) against log [AMP-PNP], where F 0 , F, and F ∞ are the fluorescence intensities of HutD alone, HutD, in the presence of various concentrations of AMP-PNP, and HutD saturated with AMP-PNP, respectively, yielded a straight line whose slope was a measure of the binding stoichiometry. The dissociation constant, K d was determined using nonlinear curve fitting analysis as per Eqs. (1) and (2). All experimental points for the binding isotherms were fitted by the least-squares method: where C 0 and C p denotes the input concentrations of the ligand and VcHutD respectively. ΔF is the change in fluorescence intensity at 334 nm (λ ex = 280 nm) for each point of titration curve and ΔF max is the same parameter when ligand is totally bound to the protein. A double-reciprocal plot of 1/ΔF against 1/(C p − C 0 ), as shown in Eq. (3) was used to determine the ΔF max .
ΔF max was calculated from the slope of the best fit line corresponding to the above plot 22 . All experimental points for the binding isotherms are fitted by least square methods using OriginPro 8.0 software.
(3) www.nature.com/scientificreports/ ATPase assay. ATPase activities of purified HutD and its mutants R18A, K44A, D166A and E167A were determined spectrophotometrically using malachite green assay by measuring the release of inorganic phosphate (Pi). As per protocol described before in Dey et al. 24 , each reaction mixture with final protein of 2.5 µM and buffer containing 50 mM Tris-HCl (pH 8.0), 300 mM NaCl, and 5 mM MgCl 2 were incubated with ATP (Sigma Aldrich), of concentration ranging from 100 to 500 μM at 298 K. Malachite Green solution was freshly prepared by adding 0.44 g of Malachite green in 0.3 M H 2 SO 4 , 2.5 ml of 7.5% ammonium molybdate, and 0.2 ml of 11% Tween 20 to a final volume of 10 ml. 200 μl of aforesaid solution was added to the reaction mixture to a final volume of 1 ml. The absorbance was measured at 630 nm within 5 min of adding the colouring reagent. Pi standard curve was prepared by using KH 2 PO 4 and plotting OD 630 against release of Pi using OriginPro 8.0 software 24 . The total Pi released by each protein upon ATP hydrolysis was obtained from the standard curve. Each protein was checked with Malachite green dye without ATP to measure the contaminant inorganic phosphate if any, and the minor absorbance thus obtained at 630 nm was subtracted from the absorbance produced by that protein upon hydrolysis of added ATP. The ATP hydrolysis without protein has served as another negative control to nullify the effect of the contaminating inorganic phosphate with ATP. All the experiments were minimally performed in triplicate.
Docking. Docking of heme with HutCD in the IF state has been performed using AutoDock Vina 32 software.
Building the lipid-protein system. The principal axis of the HutCD complex was aligned with the Z-axis, while the lipid bilayer spanned the XY plane. The HutCD complex was translated by 20 Å along the Z-axis to correctly place the TMD part within the hydrophobic region of the bilayer such that the polar head groups and hydrophilic lipid trails make contact with the specific residues of the transmembrane protein 36 The multicomponent complex was enclosed in a tetrahedral box (approximate dimensions: a = 120 Å, b = 120 Å, c = 170 Å) and solvated in explicit water with a minimum water thickness of 22.5 Å on the top and bottom of the complex. A suitable number of neutralizing ions (Na + or Cl − , depending upon the cases) were added by replacing solvent molecules using the Monte Carlo method. The number of DMPC molecules in a leaflet was in the range of 180 to 210, depending upon the protein being embedded (Fig. S3). The number of DMPC molecules on the upper and lower leaflets was slightly different (the maximum difference was 10). The area per lipid was calculated 37 by multiplying the sides of the simulation box in the XY-plane and dividing it by the average number of lipid molecules in a leaflet (Fig. S4). The coordinate of the solvated multicomponent system was saved in PDB format and was modified afterward to make it consistent with the AMBER naming scheme.
Force field for the ligands. All atom force field (FF) for the ATP cofactor was taken from the AMBER parameter database 38 maintained by the Bryce group. ATP molecules were first made fully deprotonated and then hydrogens were added using "reduce" available in the Amber Tools. All the phosphate atoms in ATP were deprotonated and the total charge was − 4e. The atom names in ATP were also changed to match the naming convention of this FF. We used antechamber 39 and MCPB.py 40 to create the FF for the heme group. In the heme group, the iron metal coordinates to four nitrogen atoms. We have used a bonded model 41 where Fe is bonded to four sp3 hybridised nitrogen. The antechamber program was used with the AM1-BCC charge method 42 to assign charge by treating heme as having a charge of − 4e and the atoms were treated with GAFF atom types 43 . The procedure we followed to treat the heme group within the AMBER software package is described with great detail in one of the Amber tutorials 41 .

MD simulation steps and parameters.
Unbiased MD simulation in this work was carried out using pmemd.cuda, the GPU version (10) of MD engine within the Amber package. The protein backbone and side chains were treated with the FF parameters from ff14SB 44 , and lipid14 FF 45 was used for the DMPC bilayer. All the simulations were carried out in explicit solvent with the three-point transferable interatomic potential (TIP3P) for waters 46 . In simulations involving only protein (i.e., the cases where bilayer was absent), a truncated octahedron box was used for solvation such that the edge of the box was at-least 10 Å away from the protein. A suitable number of Na + or Cl − ions were added to neutralize the solvated system. Periodic boundary condition was applied on all the sides of the solvated box. Bonds containing H-atoms were restrained using the SHAKE algorithm 47 which enabled us to use a time step of 2 fs for the integrator. The pressure and temperature were controlled using a weak coupling to an external bath. The solvated system was coupled to a Berendsen barostat 48 with a relaxation time of 2 ps, for maintaining a constant pressure of 1 atm. The particle mesh Ewald summation 49 was employed for calculating the electrostatic interactions part of the total Hamiltonian. We applied a distance cutoff of 10 Å for calculations of the long-range electrostatic interactions. The temperature of the system was allowed to fluctuate around a mean of 310 K by employing Langevin dynamics 50 with a collision frequency of 2/ps. Standard simulation protocols were used for the solvated protein cases 51 .
For the protein-lipid multicomponent system, the following simulation protocols were used. We followed rigorous energy minimization steps (3) to "relax" the system by eliminating any energetically unfavourable interactions that might have occurred while preparing the system using CHARMM-GUI. Minimisation was done in three steps: (1) minimised the solvent molecules holding the lipids and protein fixed by using position restraints with a force constant of 250 kcal/mol/Å 2 , (2) minimised the solvent and the lipid keeping the protein fixed and (3) finally all the atoms in the systems were set free of restraints and allowed to relax. We ran 5000 cycles of minimization steps in all three steps. The temperature of the energy minimised system was increased to 310 K linearly in two steps (0-100 K, 100-310 K). To avoid excessive solute fluctuations, we applied a weak www.nature.com/scientificreports/ restraint (10 kcal/mol/Å 2 ) on the protein and lipid residues. Langevin dynamics with a collision frequency of 1.0/ps was applied for temperature equilibration at constant volume. Each of the equilibration steps involved 100 ps of MD runs in NVT ensembles. After the heating steps, the weak restraints on the protein and lipid residues were released in steps by carrying out four MD simulations, each of 500 ps, in NPT ensembles. Monte Carlo barostat (barostat = 2) with pure semi-isotropic (ntp = 3) pressure scaling was used to achieve a constant pressure of 1 atm. To reduce buckling of the bi-layer, a constant surface tension (0.0 dyne/cm) was applied at the membrane-liquid (csurften = 3) interface 52 . Finally, we ran MD "production run" for 1000 ns with an integration time step of 2 ps in NPT ensemble.
Software used for analysis. Software used for data analysis and presentations are Origin, PyMol, UCSF chimera and Adobe photoshop.

Data availability
The pdb coordinates generated and/or analysed during the current study are not publicly available because these are generated through homology modelling. But these coordinate files will be available from the corresponding author on reasonable request. An important result of molecular dynamics simulation has been provided in the form of a movie as supplementary material. Other MD simulation results presented here will also be available in the form of movie on reasonable request.